Graphene bilayer with a twist: electronic structure 
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Electronic properties of bilayer and multilayer graphene have generally been interpreted in terms 
of AB or Bernal stacking. However, it is known that many types of stacking defects can occur 
in natural and synthetic graphite; rotation of the top layer is often seen in scanning tunneling 
microscopy (STM) studies of graphite. In this paper we consider a graphene bilayer with a relative 
small angle rotation between the layers and calculate the electronic structure near zero energy in a 
continuum approximation. Contrary to what happens in a AB stacked bilayer and in accord with 
observations in epitaxial graphene we find: (a) the low energy dispersion is linear, as in a single 
layer, but the Fermi velocity can be significantly smaller than the single layer value; (b) an external 
electric field, perpendicular to the layers, does not open an electronic gap. 



Introduction. Graphene is a two-dimensional (2D) 
carbon material, which takes the form of a planar hon- 
eycomb lattice of sp^ bonded carbon atoms. It can be 
considered as a building block for other allotropes of car- 
bon, such as graphite, fullerenes, and carbon nanotubes 
and it was first isolated by micro-mechanical cleavage of 
graphite in 2004 [H, @|- This method also produces sam- 
ples composed of two (bilayer) or more atomic layers of 
graphene (few layer graphene, FLG). FLG samples can 
also be grown epitaxially by thermal decomposition of 
the surface of SiC 

Single layer (SLG) and bilayer (BLG) graphene are 
both gapless semi-metals, if undoped, but whereas carri- 
ers in SLG have linear dispersion (leading to Dirac cones 
in energy momentum space) Q, in BLG the dispersion 
is quadratic [H|- The quantization rules for the integer 
quantum Hall effect are different for SLG 0, @, 01 and 
BLG . A controllable gap can be opened with an exter- 
nal electric field in BLG, a fact that makes it particularly 
interesting for applications (9l.[lB|. 

The properties of BLG have been interpreted under 
the assumption that the stacking of the two layers takes 
the form of AB or Bernal stacking, the most common 
in graphite. Nevertheless, AB stacking is not the only 
form of stacking found in graphite. Naturally occurring 
and synthetic (HOPG) crystals usually present a variety 
of defects which affect stacking order in the c direction. 
Turbostratic graphite is modeled by stacking graphene 
layers with random relative translations and rotations 
lll |: rotation of the top layer with respect to the bulk is 
quite common in the surface of graphite and results in the 
formation of superlattices clearly seen in STM images as 
Moire patterns jH, tl3*|. Recent detailed structural stud- 
ies of epitaxially grown FLG rule out AB stacking 
and reveal the presence of significant orientational disor- 
der of the graphene with respect to the underlying SiC 
substrate [l5| . The influence of the type of stacking on 
the electronic structure in multilayer graphene has been 



stressed by Guinea et al. |16l |. 

In this work we discuss the electronic structure of a 
bilayer with a relative, small-angle, rotation of the two 
graphene planes. We derive angles for the formation of 
periodic Moire superlattices and formulate a continuum 
electronic description in terms of massless Dirac fermions, 
coupled by a slowly varying periodic inter-layer hopping. 
We find a low energy electronic structure quite differ- 
ent from that of AB stacked bilayer, with massless Dirac 
fermions, but with a Fermi velocity {vf) substantially re- 
duced with respect to SLG. Moreover, we show that an 
external electric field does not open a gap in the electronic 
spectrum. These results are all in accord with observa- 
tions in epitaxially grown graphene, which reveal much 
the same electronic behavior as SLG in angle resolved 
photoemission ( ARPESl [ij. [itI. . transport and 
infrared spectroscopy t20| and display systematically re- 
duced values o{ vp relative to SLG |2ll |. 

Geometry. The two sublattices in layer 1 are denoted 
by A and B and in layer 2 by A' and B' . In an AB stacked 
bilayer A and B' atoms have the same horizontal posi- 
tions, iai+ja2 integers), where ai = (1/2, •\/3/2)ao, 
a2 = (— 1/2, V3/2)ao are the Bravais lattices basis vec- 
tors and oq (« 2.46 A) is the lattice constant. The SLG 
Dirac points are located at K = — K' = (47r/3, 0)/ao. 
The vertical displacement between the layers is Cq (« 
3.35 A). 

For simplicity we consider rotations of layer 2 about a 
site occupied by a B' atom (directly opposite an A atom, 
in the c direction): a commensurate structure is obtained 
if a B' atom is moved by the rotation to a position for- 
merly occupied by an atom of the same kind. The Moire 
pattern is periodic and the translation from the origin 
(center of rotation) to the B' atom's current position is 
a symmetry translation. From this we can derive a con- 
dition for the angle 9i of a commensurate rotation: 



cos 



if 



3i + 1/2 



3^2 + 3i + 1 



i = 0,1,2. 



(1) 



2 



The superlattice basis vectors are: 



ti = isLi + {i + l)a2 , 

t2 = -(?;+ l)ai + (2i + l)a2, 



(2) 



(i = is an AA stacked bilayer). 
of the superlattice is L = |ti| = v3i^ 



The lattice constant 
3i + lao. STM 
measurements of the surface of graphite observed 
superlatices with periods of L = 66 A and angles = 2.1° 
corresponding to j = 15 in ([1]) and a unit cell with 2884 
atoms, making ab initio descriptions rather impractical. 
The reciprocal lattice vectors are: 



Gi = 



47r 



[(3i + l)ai +a2] 



(3) 



3 (3i2 + 3i + l) 
G2 - 3(3,.'l3, + i) K3. + 2)a, + (3^ + l)a2].(4) 

Continuum description. The Hamiltonian for the bi- 
layer with a twist has the form Hi +Ti.2+ 'H± , with the 
intra-layer Hamiltonian, H1+H2, given by (we use units 
such that h—1): 



Ui = ~t'^c\(v.i)[cB{Yi 



So) 



CB(rj + So - ai) 



CBi^i + So - a2)] + h.c. 



(5) 
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+ CA'{rj 



h.c. 



(6) 



where Ca{v) is the destruction operator for the state in 
sublattice a at horizontal position v; a = A, B \n layer 1 
and a = A' , B' in layer 2; a.[ and a.2' are obtained from ai 
and a2 by a rotation by 6 about the origin; = ■ma.i+na.2 
for Hi and rj = ra'^ + sa'2 for H2 {m,n,r,s, integers); 
So = (ai + a2)/3 and Sq = (a'^ + a2)/3. 

To study the low energy spectrum near the K (K') 
point, we go to the continuum limit, with the standard re- 

1 /2 

placement Ca(r) v/ -01, 0(1") exp(iK • r) where ipi.a{r) 
is a slowly varying field on scale of the lattice constant 
{vc is the unit cell volume). Due to the rotation, the wave 
vector in layer 2 is shifted to = 47r(cos 0, sin0)/(3ao), 
so Ca'(r) Dy^'!/'2,a(r) exp(iK^ • r). For small angles of 
rotation the modulation of inter-layer hopping has a long 
wavelength and the coupling between different valleys (K 
and K') can be ignored. Hence, in the long-wavelength 
limit the decoupled Hamiltonian can be written as: 



k 

"^2 = V'2,k'^T^-k7/'2,k, 



(7) 
(8) 



where vp = and t = {Tx,Ty) are Pauli matrices. 

The coordinate axes have been chosen to coincide with 
those of layer 1, so the Hamiltonian of layer 2 involves a 



extra rotation by 9, the angle between the two layers and 
T« =e+*«"=/2(T„Ty)e-'^^^/2. 

To model the inter-layer coupling, H±, we retain hop- 
ping from each site in layer 1 to the closest sites of layer 
2 in either sub-lattice. We denote by "(r) the hori- 
zontal (in-plane) displacement from an atom of layer 1, 
sub-lattice a{a — A, B) and position r, to the closest 
atom in layer 2, sub-lattice f3' {f3' — A',B'). Denoting 
by t± (6) the hopping between pz orbitals with a relative 
displacement co -I- S, one gets 

^-L = E 4(r^)c/3' (r, +<5^'"(r,)) + h.c. 

(9) 

where t± ((5"''(r)) = t'^{r), is the inter-layer, position 
dependent, hopping between orbitals with a relative 
displacement Co -I- S; AK = — K is the relative shift 
between corresponding Dirac wavectors in the two layers; 
4>i,k,a = ^i,k±AK/2,a IS the Fourier component of -tpi^a (r) 
for momentum k± AK/2, the plus sign applying in layer 
1 and the minus one in layer 2. With this choice, the 
Dirac fields (j)i,k,a with the same k vector in both lay- 
ers correspond to the same plane waves in the original 
lattice; the Dirac cones occur at k = —AK/2 in layer 1 
and AK/2 in layer 2. Replacing the operators in eq. ((9| 
with the Dirac fields, using ipi^pir + S^°'{r)) « '0i^^(r), 
since the Dirac fields are slowly varying on the lattice 
scale, and Fourier transforming, the low energy effective 
Hamiltonian, near K, is 

k,al3 ^ ^ 

+ VfY, 'l'ik,a'^% ■ ( k - ^ ) 02,fc,/3 
k,a,l3 ^ 
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a, 13 k,G 



if (G) 



^lfc+G,a'^2A/3+h.C. . (10) 



For commensurate structures, the function 
t"'^(r) exp (iK^ • (5"^(r)) is periodic and has nonzero 
Fourier components only at the vectors G of the 
reciprocal lattice : 



tf(G) 



1 

Vc 



■S°"^{r)^-iG-r 



(11) 



The integral is over the unit cell of the superlattice and it 
will ultimately be calculated by a sum over the sites of the 
Wigner-Seitz unit cell since t"*^ (r) is only defined at those 
points. This Hamiltonian describes two sets of relativistic 
Dirac fermions (with shifted degeneracy points) coupled 
by a periodic perturbation. 

Fourier amplitudes of inter-layer coupling. To deter- 
mine the hopping t±{5) as a function of the horizontal 
shift S"'^{r) we express it in the Slater- Koster parameters, 
Vppa{d) and VppTr{d), where d is the distance between the 
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Figure 1: First Brillouin zone (FBZ) of the super-lattice cen- 
tered at mid-point between Dirac points K and K*. Note 
that the zero energy states of the two layers, k = — AK/2 
and k — AK/2, marked with ®, are half-way to the zone 
boundary; AK is a vertex of the FBZ. 
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Table I: The most important Fourier amplitudes are shown 
in this table (all others are smaller by at least a factor of 
5 for angles smaller than 10°). The first and second lines 
express exact results: t± is real. In the last two lines there 
are corrections to these results of order ao /L where L is the 
period of the superlattice. 



■\/cq + (5^ . For the d dependence 



pr>a{d) and VppTt{d) we used the parametrization of 



two atomic centers, d 
oiVp 

ref. 14'p7r(ao/%/3), is the in-plane nearest neighbor 

hopping, t, and Vppcrico) if the inter-layer hopping, t±, in 
an AB stacked bilayer. The contribution of T^p^ turns 
out to be negligible and t\^{5) is proportional to t^: for 
5 = ao/%/3, t^{5)/tx_ w 0.4. We have calculated (5"'3(r) 
numerically for any angle of rotation. Using various sym- 
metries and relations valid in the limit L (small 
angles) we were able to derive the results of Table [H The 
values of i^'^{G) are equal and real, by symmetry, for 
G = 0, G = — Gi and G = — Gi — G2 and much smaller 
for all other G vectors. The remaining Fourier ampli- 
tudes can be expressed in terms of tf^{G). 

Results and discussion. In the absence of the inter- 
layer coupling, 1-L±, states with energy close to zero occur 
at k = -AK/2 in layer 1 and k = -l-AK/2 in layer 2. 
The results of Table U imply that the states of momen- 
tum k in layer 1 are coupled directly only to states of 
layer 2 of momentum k, k 4- Gi and k + Gi 4- G2; con- 
versely the states of momentum k in layer 2 only couple 
to states k, k — Gi and k — Gi — G2. To investigate 
the spectrum at a momentum k close to zero energy, 
we truncated the Hamiltonian to include only these six 
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Figure 2: The energy of the two states with smaller \ek \ for 
6 = 3.9° (« = 8); panel (a): k varying form - AK to AK (two 
vertexes of the FBZ) along the line passing the degeneracy 
points, —AK/2 to AK/2; panel (b): along a line parallel to 
G2 passing AK/2. 



momentum values (three for each layer) giving a 12 x 12 
matrix to diagonalize. The geometry of the first Brillouin 
zone (FBZ) of the superlattice (fig. [l]) implies that the 
states near the degeneracy point in either layer couple 
only to states of energies -iivp^K — zLvpK x 2sin(0/2) 
where AK = |AK| and K = 47r/(3ao). This turns out 
to be the essential difference between this problem and 
that of the unrotated bilayer. In the latter, the degener- 
acy points of both layers occur at the same momentum 
and the inter-layer hopping couples two doublets of zero 
energy states. In the present case we have one doublet of 
zero energy states coupling to three pairs of states at fi- 
nite energies, zLvpAK. As a result, the linear dispersion 
near zero energy is retained. In Fig. [2] we plot the ener- 
gies of the states with smallest \ek\ alongtwo lines in the 
FBZ; the parameters are ij_ = 0.27 eV^ and = 3.9° 
(i = 8, L = 36 A), which give vpAK « 0.76 eV and 
i± = 0.11 eV. 

The persistence of the Dirac cones can be understood 
by considering the limit where t±/ (vpAK) ^ 1 (in the 
situation represented in fig. El i±/{vFAK) « 0.14). Con- 
sider, for instance, the vicinity of the degeneracy point 
of layer 1, k = —AK/2 -I- q. It is clear that the Hamil- 
tonian H(k) has the form H(k) = H(-K/2) -h V(q) with 
V(q) linear in q. In H(— K/2), which contains the inter- 
layer coupling, the doublet at zero energy couples with 
an amplitude ~ i± to six states (of layer 2) with ener- 
gies ±vfAK. Using perturbation theory one can derive 
an effective Hamiltonian in the space of the zero energy 
doublet by considering the mixing of these six states in 
layer 2 to first order in i±/ {vpAK). The degeneracy 
is not lifted, although there is a small shift in energy, 
eo = &t\sa\{9/2)/{vFAK). For small q we can treat 
V(q) as a perturbation in the subspace of this doublet: 
the effective Hamiltonian matrix has the form character- 
istic of a Dirac cone 



He 



eo vpq 
vpq eo 



with q = qx + iqy To second order in t^^/vpAK, 
the renormalized Fermi velocity is given hy vp/vp = 



4 




-1 -0.5 0.5 1 



k/A K 

Figure 3: The energy of the two sates with smaller |efe| in 
the presence of a potential difference V = 0.3 V between layers; 
k varies from — AK to AK (two vertexes of the FBZ) along 
the line passing the degeneracy points, — AK/2 to AK/2; the 
remaining parameters are the ones used in Fig. [21 

1 — 9 (t±/{vF'!^K))'^ . This significant depression of 
the value of the Fermi velocity vf relative to the value 
of single layer graphene is a tell-tale sign of the pres- 
ence of a bilayer with a twist. The perturbative results 
slightly overestimates the downward renormalization of 
vp, because of the contributions of higher order terms in 
i±/ (vf^K), especially at smaller angles (smaller AK, 
larger i±/ {vpAK)). Ref. [2l| reports several observa- 
tions of values ofvp in the range 0.7 ^ 0.8 x lO^ms"^, in 
epitaxial graphene, 20 to 30% lower than in single layers. 

Another important consequence of the rotation be- 
tween layers occurs when there is an electric poten- 
tial difference between layers. To the Hamiltonian ifTO]) 
this adds a term Vext = -(^/2) Sfc.a </'i,/c,a'Ai,fc," + 
{V/'^)J2k,a'l'2,k,a^'^,k,a- It is kuowu that in the AB 
stacked bilayer a gap opens in the spectrum in the pres- 
ence of an external electric field between layers [1, [lB |. 
However, it is clear from the discussion above that the 
cones present in the bilayer with a twist are essentially 
the Dirac cones of each layer perturbed by the admix- 
ture of states of the opposing layer, which are distant in 
energy. As such, we expect that a potential difference 
between the layers, V, should merely give rise to a rela- 
tive shift of the energies of the degeneracy points in each 
cone, at least as long asV< vpAK. This expectation is 
borne by the results shown in Fig. [Sj the Dirac cones are 
shifted but there is no gap in the spectrum. 

The results of Table U imply that a small angle rotation 
destroys the particle-hole symmetry of an AB stacked bi- 
layer (with hopping only between A and B' atoms) . The 
Fermi level of an undoped sample need no longer be at 
zero energy; in turbostratic graphite, for instance, it is 
shifted to 0.11 eV jUj. This calculation, being Hmited 
to energies close to zero, cannot determine the absolute 
position of the Fermi Level as a function of carrier con- 
centration. 

In conclusion, we presented a detailed geometrical de- 



scription of a bilayer with a relative rotation between the 
layers. We developed a continuum description valid for 
small angles of rotation and analyzed the energy spec- 
trum close to zero energy. We found that the Dirac 
cones of a single layer graphene remain present in the 
bilayer, but with a significant reduction of the Fermi ve- 
locity especially for very small angles of rotation. A new 
energy scale is introduced vpAK — vpK x 2sin(6'/2) 
where K — 47r/(3ao) and 6 is the angle of rotation; the 
dispersion relation is only linear for energies such that 
|efc| < vpAK. Unlike the case of the AB stacked bilayer, 
a potential difference between layers does not open a gap 
in the spectrum. These results show that a small stacking 
defect such as a rotation can have a profound effect on 
the low energy properties of the bilayer and are in accord 
with several observations in epitaxial graphene. 
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